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Abstract 

Pressure (density) and velocity boundary conditions inside a flow domain are studied for 
2-D and 3-D lattice Boltzmann BGK models (LBGK) and a new method to specify these 
conditions are proposed. These conditions are constructed in consistency of the wall boundary 
condition based on an idea of bounceback of non-equilibrium distribution. When these con- 
ditions are used together with the improved incompressible LBGK model Q, the simulation 
results recover the analytical solution of the plane Poiseuille flow driven by pressure (density) 
difference with machine accuracy. Since the half-way wall bounceback boundary condition is 
very easy to implement and was shown theoretically to give second-order accuracy for 2-D 
Poiseuille flow with forcing, it is used with pressure (density) inlet/outlet conditions proposed 
in this paper and in Q to study the 2-D Poiseuille flow and the 3-D square duct flow. The 
numerical results are approximately second-order accurate. The magnitude of the error of 
the half-way wall bounceback is comparable with that using some other published boundary 
conditions. Besides, the bounceback condition has a much better stability behavior than that 
of other boundary conditions. 

1 Introduction 

The lattice Boltzmann equation (LBE) method has achieved great success for simulation of trans- 
port phenomena in recent years. Among different LBE methods, the lattice Boltzmann BGK 
model is considered more robust |J. Some recent theoretical discussions on LBGK ||, ||] have 
enhanced our understanding of the method and the effect of boundary conditions. They explains 
why the velocity boundary condition for the 2-D triangular LBGK model proposed in || generates 
results of machine accuracy for plane Poiseuille flow with forcing, and the bounceback or equilib- 
rium scheme generate an first-order error to the velocity. Moreover, the bounceback scheme with 
the wall located half-way between a flow node and a bounceback node (it will be called "half-way 
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wall bounceback" thereafter) is shown theoretically to produce results of second-order accuracy 
for the simple flows considered. 

The above mentioned results with Poiseuille flows are obtained with external forcing to drive 
the flow. In practice, however, a flow is often driven by pressure difference, and the pressure 
gradient in many cases cannot be replaced by an external force in LBGK computations. In this 
situation, boundary conditions usually need be implemented by giving prescribed pressure or 
velocity on some "flow boundaries", which are not solid walls or interfaces of two distinct fluids. 
Instead, they are imaginary boundaries inside a flow domain (e.g. inlet and outlet in a pipe flow). 
Their existence is purely for the convenience of study. The implementation of these boundary 
conditions in LBGK is very important but it has not yet been well studied. 

In lattice Boltzmann method, a specification of pressure difference amounts to a specification 
of density difference. Early works (see, for example, @) to implement pressure (density) flow 
boundary condition is simply to assign the equilibrium distribution computed with the specified 
density and some velocity (maybe zero) to the distribution function. This method introduces 
significant errors. Skordos || proposed to add a term to the equilibrium distribution to improve 
it, the scheme requires the gradient of density and velocities at boundaries. Inamuro et al. § and 
Maier et al. [ |l0| proposed new boundary conditions for lattice Boltzmann simulations. In their 
simulation of Poiseuille flow with pressure (density) gradient, the pressure boundary condition 
was treated in different way compared to their wall boundary condition. Chen et al. 0] also 
proposed a way to specify general boundary conditions including flow boundary conditions. 

In this paper, we propose a way to specify pressure or velocity on flow boundaries. They are 
treated as one type of general boundary conditions, which are based on an idea of bounceback 
of non-equilibrium part with modifications. When applied to the modified LBGK model, These 
boundary conditions produce results of machine accuracy for 2-D Poiseuille flow with pressure 
(density) or velocity inlet/outlet conditions. It is also noticed that although all the proposed new 
boundary conditions (||, |(| ||, [R]]) including the boundary conditions in this paper yield improved 
accuracy compared to the bounceback boundary condition, they are difficult to implement for 
general geometries, because there is a need to distinguish distribution functions according to their 
orientation to the wall. Besides, there are additional works or different treatments at corner 
nodes. On the other hand, the complete bounceback scheme does not distinguish distribution 
functions and is very easy to implement in a parallel way, which was considered as one of the 
advantages of LGA or LBE method. In this paper, the half-way wall bounceback boundary 
condition with two flow boundary conditions are applied to the 2-D Poiseuille flow and a 3-D 
square duct flow using the d2q9i and d3ql5i lattice Boltzmann models respectively. The results 
are approximately second-order accurate. The error is comparable with that using some published 
boundary conditions. Moreover, the half-way wall bounceback boundary condition for stationary 
walls is much more stable than the boundary conditions in || ^, ||, [l(J and in this paper. Thus, we 
recommend to use the half-way wall bounceback boundary condition for stationary walls and use 
the boundary conditions proposed in this paper or in [g, only for flow boundary conditions. 



2 Governing Equation 

The square lattice LBGK model (d2q9) is expressed as (Q,||,||]): 

/i(x + 5e u t + 5) - /i(x, t) = -i[/i(x, t) - /f 9) (x, t)}, i = 0, 1, 8, 



(1) 
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where the equation is written in physical units. Both the time step and the lattice spacing have 
the value of 8 in physical units. /j(x, t) is the density distribution function along the direction ej 
at (x,i). The particle speed e^'s are given by ej = (cos(-7r(i — l)/2), sin(7r(z — l)/2),i = 1,2,3,4, 
and ej = \/2(cos(7r(z — 4 — i)/2), sin(-7r(i — 4 — i)/2),z = 5,6,7,8. Rest particles of type with 
eg = is also allowed (see Fig. 1). The right hand side represents the collision term and r is the 
single relaxation time which controls the rate of approach to equilibrium. The density per node, 
p, and the macroscopic flow velocity, u = (u x ,u y ), are defined in terms of the particle distribution 
function by 

8 8 

fi = p> J2 fo ei = p u - ( 2 ) 

4=0 4=1 

The equilibrium distribution functions f^ eq \x,t) depend only on local density and velocity and 
they can be chosen in the following form (the model d2q9 [12]): 

if 9) =i i p[l + 3(e i -u) + ~(e i -u) 2 -~u-u], t = -, ^ = 1^ = 1:4; i i = -L i = 5:8. (3) 

A Chapman-Enskog procedure can be applied to Eq. (jl|) to derive the macroscopic equations of 
the model. They are given by: the continuity equation (with an error term 0(5 2 ) being omitted): 

^ + V-(pu)=0, (4) 

and the momentum equation (with terms of 0(5 2 ) and 0(5u 3 ) being omitted): 

dt(pu a ) + d j3 (pu a up) = -d a (c 2 s p) + dp(2vpS aj 3), (5) 

where the Einstein summation convention is used. S a p = \{d a Uf3 + dpu a ) is the strain-rate tensor. 

2 1 2r - 1 

The pressure is given by p = c^p, where c s is the speed of sound with c s = -, and v = o, 

3 6 

with v being the the kinematic viscosity. The form of the error terms and the derivation of these 
equations can be found in 0, 15]. 

For 2-D case, we will take the Poiseuille flow as an example to study the pressure (density) or 
velocity inlet/outlet condition. The analytical solution of Poiseuille flow in a channel with width 
2L for the Navier-Stokes equation is given by: 

y 2 dp dp 

u x = u {l-j^), u y = 0, ^ = "G, qZ = ^ ( 6 ) 

where the pressure gradient G is a constant related to the centerline velocity u$ by 

G = 2puu /L 2 , (7) 

and the flow density p is a constant. The Reynolds number is defined as Re = uq{2L)/v. 

The Poiseuille flow is an exact solution of the steady-state incompressible Navier-Stokes equa- 
tions with constant density po: 

V • u = 0. (8) 

p 

dp(u a U{3) = -d a (—) + vdppua, (9) 
Po 
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On the other hand, the steady-state macroscopic equations of the LBGK model are different from 
the incompressible Navier-Stokes equations Eqs. (|§||^) by terms containing the spatial derivative 
of p. These discrepancies are called compressibility error in LBE model. Thus, when pressure 
(density) gradient drives the flow, u x in a LBGK simulation increases in the x— direction, the 
velocity profile from the simulation is no loger a parabolic profile. For a fixed Mach number (uq, 
fixed), as 5 — > 0, the velocity of the LBGK simulation will not converge to the velocity in Eq. @ 
because the compressibility error becomes dominant. This makes the comparison of u x with the 
analytical velocity of Poiseuille flow somehow ambiguous. 

To make a more accurate study for Poiseuille flow with pressure (density) or velocity flow 
boundary condition, it is better to use the improved incompressible LBGK model proposed in 
M. The model (called d2q9i) is given by Eq. @ with the same e, and the following equilibrium 
distributions: 

/f ?) =t 1 [/ ) + 3e 1 -v + ^(e ! .v) 2 Jvv], t Q = -, = = 1 : 4; ti = ^,i = 5:8. (10) 
and 

i=0 i=0 i=l i=X 

where v = (v x ,v y ) (like the momentum in the ordinary LBGK model) is used to represent the 
flow velocity. The macroscopic equations of d2q9i in the steady-state case (apart from error terms 
ofO(5 2 ): 

V-v = 0, (12) 

dp{v a vp) = -d a (c 2 s p) + vdpf3V a , (13) 

are exactly the steady-state incompressible Navier-Stokes equation with constant density po- In 
this model d2q9i, pressure is related to the calculated density by c 2 p = p/po (c 2 = 1/3), and 
v = 2r r 5. The quantity p/po will be called the effective pressure. Although the macroscopic 
equations of d2q9i in the steady-state case has an error of 0(5 2 ) to the steady-state Navier-Stokes 
equation, for some special flows like the Poiseuille flow, it is possible that this error disappears 
with suitable boundary conditions. 



3 Pressure or Velocity Flow Boundary Condition of the 2-D 
Square Lattice LBGK Model 

In this section a new boundary condition is proposed based on an idea of bounceback on non- 
equilibrium part as follows: take the case of a bottom node in Fig. 1, the boundary is aligned with 
x— direction with f±, fy, /g pointing into the wall. After streaming, fo, /i, /3, /4, /V, fs are known. 
Suppose that u x , u y are specified on the wall, we want to use Eqs. @ to determine /2, /s, fe and 
p (originated in H), which can be put into the form: 

h + h + k = p - Ob + fx + h + U + h + fs), (14) 

h-k = pu x - (fi -h-h + fs), (15) 

f2 + h + k=pUy + (/4 + fl + fs)- (16) 
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Consistency of Eqs. ( |l4| , |lq ) gives 

P = r^—ifo + fa + fa + 2(/ 4 + fa + fa)]- (17) 

1 — Uy 

We assume the bounceback rule is still correct for the non-equilibrium part of the particle 
distribution normal to the boundary (in this case, fa — = fa — f^)- With fa known, fa, fa 
can be found, thus 

fa = fa + T^PUy, 

fa = fa-^(fl- fa) + ^PUx + ^PUy, 
1 11 

fa = fa + 2 (A ~ /s) - 2 pUx + Q pu y ( 18 ) 

The collision step is applied to the boundary nodes also. For non-slip boundaries, this bound- 
ary condition is reduced to that in juj. A detailed discussion of implementation of boundary 
conditions on stationary walls in 3D case was given in [1C]. 



3.1 Specification of Pressure on a Flow Boundary 

Now let us turn to the pressure (density) flow boundary condition. In M, 10], the pressure (density) 



boundary condition was treated in different way compared to the wall boundary condition. Chen 
et al. use an extrapolation scheme on an additional layer beyond a boundary to determine 
the incoming fas before the streaming step, their treatment of pressure boundary condition is 
consistent with the wall boundary condition. In this paper, we treat pressure (density) flow 
boundary condition the same way as the velocity wall boundary condition. Its derivation is based 
on Eq. (0) as for velocity wall boundary condition. Suppose a flow boundary (take the inlet in 
Fig. 1 as example) is along the y— direction, and the pressure (density) is to be specified on it. 
Suppose that u y is also specified (e.g. u y = at the inlet in a channel flow). After streaming, 
fa, fa, fa, fa, fa are known, p = pi n ,u y = are specified at inlet. We need to determine u x and 
fa, fa, fa from Eq. (§) as following: 

fl + fa + fa = Pin ~ (fa + fa + fa + fa + fa + fa), (19) 

fa + fa + fa = p in u x + (/ 3 + fa + fa) , (20) 
fa-fa = -fa + fa~fa + fa- (21) 
Consistency of Eqs. ( |l9| , pG| ) gives 

Ux = 1 - [/o + A+A+ 2 (A + & + (22) 

Pin 

We use bounceback rule for the non-equilibrium part of the particle distribution normal to 
the inlet to find fa — f[ eq ^ = fa — ' ■ With fa known, fa, fa are obtained by the remaining two 
equations: 

fa = fa+ T^PinUx, 

fa = fa ~ \{fa ~ fa) + \pinU x , 

1 1 

fa = fa+ ^(fa - fa) + -PinUx- (23) 

2 o 
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The corner node at inlet needs some special treatment. Take the bottom node at inlet as 
an example, after streaming, /3,/4,,/V are known; p is specified, and u x = u y = 0. We need to 
determine fi, f'2, f'5, k, h- We use bounceback rule for the non-equilibrium part of the particle 
distribution normal to the inlet and the boundary to find: 

h = h + U[ eq) - ft ] ) = h, h = h + (ft ] - A eq) ) = h, (24) 

Using these /i,/a in Eqs. (|20|J?ll) , we find: 

h = /7, k = k = \[Pin - (/o + fl + h + h + h + k + h)\ (25) 
Similar procedure can be applied to top inlet node and outlet nodes including outlet corner nodes. 

3.2 Specification of Velocity on a Flow Boundary 

In some calculations, velocities u x ,u y are specified at a flow boundary (take the inlet in Fig. 1 as 
example). In the flow region of the inlet or outlet, this is actually equivalent to a velocity wall 
boundary condition and can be handled in the same way as given at the beginning of the section. 
The effect of specifying velocity at inlet is similar to specifying pressure (density) at inlet. Density 
difference in the flow can be automatically generated by the velocity inlet condition. 

At the inlet bottom (non-slip boundary), special treatment is needed. After streaming, 
fli f%-> k, ki /s need to be determined. Using bounceback on normal distributions gives: 

fl = f3, /2 = fi- 

Expressions of x, y momenta give: 

h — /e + h = -(h - h - h) = 

h + k-h = -{h-h-ft) = h, (26) 

or 

/5 = fit 

k = k = \[p-{k + h + h + k + h + k + h)l (27) 

but there is no more equation available to determine p. The situation is similar to a corner wall 
node (the intersection of two perpendicular walls). In this situation, since p is expected to be 
constant at the inlet, p at the inlet bottom node can be taken as the p of its neighboring flow 
node, thus the velocity inlet condition is specified. It is noted that this treatment is only for 
a special study to produce the analytical solution with model d2q9i. In practice, however, the 
half-way wall bounceback is recommeded as boundary conditions for stationary walls and there 
is no need for special treatment at corner nodes. This will be discussed in section 4. 

From the discussion given above, we can unify boundary conditions (on a wall boundary or in 
a flow boundary) in 2-D simulation on a straight boundary as: 

• Given u x ,u y , find p and unknown /j's. 
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• Given p and the velocity along the boundary, find the velocity normal to the boundary and 
unknown /j's. 

The above discussion is for flat flow boundaries aligned with a plane spanned by two particle 
velocities of type I. It is not the purpose of this paper to derive a flow boundary condition in 
a general geometry. If a flow boundary has a complicated geometry or if it is not aligned with 
lattice directions, schemes based on extrapolations like the ones in |2j, can be used. However, 
as pointed out in on a convex edge or at a convex corner, there are too few unknown /j's, if 
the available /j's are used, then any choice of the unknown /j's may not make the velocity correct. 
There have not been analytical or numerical studies on the order of accuracy for these situations. 
Thus the order of accuracy is not clear for these cases. Hence, the equilibrium scheme can be 
considered as well. The equilibrium scheme at boundaries gives second-order accuracy if r = 1 



but only first-order when r / 1 [g, 16 1. 



3.3 Boundary Conditions for the Modified Incompressible Model d2q9i 

The velocity wall boundary condition and flow boundary conditions for d2q9i are similar to that 
of d2q9. The derivation is based on equations Yh=o fi = P an d Z)f=i e ifi = v an d hence some 
modifications are needed as follows: 



In wall boundary condition, Eq. ( |17| ) is replaced by 

P = v y + [/o + fx + h + 2(/ 4 + h + h)\ (28) 



and in Eq. (|18[) , pu x ,pu y are replaced by v x ,v y respectively. 
In pressure flow boundary condition, Eq. (22) is replaced by 



[/0 + /2 + /4 + 2(/ 3 + /6 + / 7 )], (29) 



and in Eq. fl23|), pi n u x is replaced by v x . 



3.4 Numerical Results of Model d2q9i 

We report and discuss the numerical results for Poiseuille flow with our wall and flow boundary 
conditions. The simulation is performed on the model d2q9i. The main result in the simulation 
is the achievement of machine accuracy. The width of the channel is assumed to be 2L = 2. We 
use nx,ny lattice nodes on the x— and y— directions, thus, 5 = 2/(ny — 1). The initial condition 
is to assign /j = f^ eg ' computed using a constant density po, and zero velocities. The steady-state 
is reached if 

EiE 7 - \v x {i,j,t + 5) - v x (i,j,t)\ + \v y (i,j,t + 5) - v y (i,j,t)\ 



< 5 ■ Tol. (30) 



Tol is a tolerance set to 10 12 in this section. 

We also define a maximum relative error of velocity (v x ,v y ) as in [16] 



err m = max J , (31) 

u 
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where u^Uy is the analytical velocity, and uq is the peak velocity, the maximum is throughout 
the flow. 

For model d2q9i, we carried out simulations with a variety of Re, nx, ny, uo using the pressure 
or velocity flow boundary condition. All simulations in this paper use double-precision. The 
range of Re is from 0.0001 to 30.0; the range of r is from 0.56 to 20.0 and the range of uq is from 
0.001 to 0.4; the largest density difference simulated (not the limit) is pi n = 5.6, p ou t = 4.4 with 
nx = 5, ny = 3 corresponding to an effective pressure gradient of G' = 0.1, where G' is defined as 
G' = -ijg. The magnitude of average density po is 5, but it is irrelevant for the simulation [jj]]. 

For all cases where the simulation is stable, the steady-state velocity and density show: 

• The velocity field v x is uniform in the x-direction, it is accurate up to machine accuracy 
compared to the analytical solution in Eq. v y is very small with maximum of \v y \ in the 
whole region being in the order of 10 -13 . For example, for nx = 5,ny = 3, uq = 0.1, r = 0.56, 
Re =10, the maximum relative error of velocity is 0.1816- 10 -11 , while the maximum relative 
error of density is 0.3553 • 10~ 15 , and the maximum magnitude of v y is 0.5551 • 10~ 15 . The 
results for other cases are similar to this example. 

• The density is uniform in the cross channel direction, and linear in the flow direction, the 
computed value and the analytical value of the density gradient differ only at the 14th digit. 

It is also noticed that with pressure (density) gradient to drive Poiseuille flow, the maximum 
Reynolds number which makes the simulation stable is far less than that with external forcing, 
which is also reported in (2|. 

Similar results with machine accuracy are obtained by specifying the analytical velocity profile 
given in Eq. ([]) at inlet and pressure (density) at outlet by using the flow boundary conditions in 
this paper. In the case, there is a uniform pressure (density) difference in the region. The value 
of the density difference depends on uq and the outlet density. 

The accurate results in the model d2q9i give us confidence about the flow boundary conditions 
proposed. 



4 On Half-way Wall Bounceback Boundary Condition 

4.1 Reconsider Bounceback 

The bounceback rule with wall placed at the bounceback nodes gives an first-order error to 
velocity, both at the boundary and throughout the flow. This has been showed analytically || 
for some simple flows and computationally [15| for 2D cavity flows. The same can be said for 
the flows studied in this paper based on numerical results. There have been efforts to replace 
the bounceback boundary condition with new boundary conditions for LBE simulation. Various 
boundary conditions || |(| [8|, [9|, [l(], 16 1 including the one in this paper are proposed. While these 



boundary conditions indeed improve the accuracy for some simple flows over the bounceback 
scheme, extra works are still needed to apply them to a general situation. First, for example, at 
a concave edge [1C], the boundary conditions in || || || [R], [Uj] and in this paper do not give 



a natural way to determine the density, and since there are more unknowns at a node than in 
a plane wall node, some additional assumptions are needed to apply the boundary condition. 
Second, all the boundary conditions in |2], ||, ||, Q [l(], 16] including the one in this paper require 



different treatments on /j's depending on their orientation to the wall, the implementation of 
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these boundary conditions is formidable for complicated geometry like that in a porous media or 
for a situation where a wall is not aligned with any of the /j's direction. On the other hand, the 
"complete" bounceback scheme which assigns each /j the value of the fj of its opposite direction 
with no relaxation on the bounceback nodes is very easy and convenient to apply, the treatment 
is independent on the direction of /j's, which is one of the major advantages of LGA or LBE 
method. Several theoretical studies have shown that if the wall is placed at the half way between 
the bounceback row and the first flow row ("half-way wall bounceback"), the scheme gives a 
second-order accuracy @, |, 0, Hi for some simple flows including an inclined channel flow and a 
plane stagnation flow. For example, if the d2q9 model with the half-way wall bounceback is used 
to simulate the 2-D Poiseuille flow with forcing, the error of the velocity (it is the same for any 
node) is given by || 

^-^ = -^M4r-5 ) + 3]^ 2 (32) 
3 

where u^Uj is the analytical and computed x- velocity respectively, uq is the center velocity. For 
a fixed r, the error is second-order in the lattice spacing 5. Of course, large value of r will give 
large errors, which are reported in some papers H M, giving an impressed vision how bad the 
bounceback boundary condition was. However, for practical purposes, there is no need to take 
large value of r in a simulation. Besides, the Chapman-Enskog procedure does not work right 
when r is large. For r between 0.5 and 1.25, the magnitude of the error given in ( p2| ) is less than 
or equal to 1.1 uq5 2 . Thus, it is worthwhile to consider this boundary condition in more general 
situations especially in 3-D flows and some results are reported in section 4.2, 4.3 and section 5. 



4.2 Results of Model d2q9i 

To reduce the effect of compressibility error, the model d2q9i is used. Simulations with model 
d2q9 were also carried out, the results have very similar behavior for order of convergence but the 
magnitude of error in d2q9 is greater. We use d2q9i to simulate Poiseuille flow with pressure or 
velocity flow boundary condition. For the half-way wall bounceback, if there are ny nodes on the 
y— direction, then the first and last nodes are the bounceback nodes with the wall being located 
half-way between the bounceback node and the first flow node, there are ly = ny — 2 lattice steps 
across the channel and the lattice spacing is 5 = 2/(ny — 2) = 2/ly (in the case of the boundary 
condition in |@] 8 = 2/(ny — 1) = 2/ly). The length of the channel is set to twice as the width. 
At the inlet/outlet, the bounceback is also used at the nodes on the bounceback rows, thus, there 
is no additional treatment for the corner nodes at the inlet/outlet. 

Two inlet/outlet (I/O) conditions with the half-way wall bounceback were tested: 

1. I/O No. 1: the flow boundary condition given in this paper. 

2. I/O No. 2: the flow boundary condition proposed in M, It assumes an additional layer 
of nodes beyond the boundary flow nodes and uses an extrapolation formula to derive the 
incoming /j's of the additional layer before streaming. 

The I/O and boundary condition in [Q] were also used to compare result with the the half-way 
wall bounceback. 

For the study of accuracy, we fix r and the Reynolds number (which affects stability) and 
halve 5 each time and calculated the error in the result. Simulations with different Re and 
r < 1.3 are performed. Three examples with pressure specified at inlet and outlet are reported 
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on Table I, The example uses three sets of parameters: (1) r = 0.6, Re = 10 , (2) r = 0.8, Re 
= 10, and (3) r = 1.1, Re = 1 (Re is restricted so that uq for the smallest number of lattice 
nodes is still small), and po = 5. The quantity Tol in Eq. (|3C|) is set to 10~ 8 . As 5 halves, uq 
also halves, so does the Mach number, the same way as in the study of duct flow in flO|| . As 5 
changes, v and the effective pressure gradient G' and then pressure (density) at inlet/outlet also 
changes. Ix, ly are used to represent the number of lattice steps in x— and y— directions, we use 
Ix = 8, 16, 32, 64, 128, ly = 4, 8, 16, 32, 64 respectively to do the simulation. 

The convergence result is summerized in Table I. The ratio of two consecutive maximum 
relative errors is also shown. The order of convergence from a least-square fitting (a linear least- 
squares fitting to logarithms of error and 5) is shown in the last column. 

For the cases of I/O Nos. 1, 2 with the half-way wall bounceback, the maximum relative 
velocity errors are very close, the ratio is approximately equal to 4, indicating a second-order 
accuracy. The magnitude of errors in the half-way wall bounceback is close to (sometimes smaller 
than) the result using I/O and boundary condition in (2j. For the cases of the half-way wall 
bounceback, it is also observed that 

• Velocity v x is generally uniform in the x— direction. 

• v y is small compared to uq, with maximum of \v y \ being less than 0.011-uo in all cases. As ly 
increases, the maximum of \v y \/uo decreases. I/O No. 2 gives much smaller max|f y | than 
I/O No. 1. 

• The density is approximately uniform in the cross channel direction, and linear in the flow 
direction. The maximum relative density error is much less than the maximum relative 
velocity error and it decreases much faster than the maximum relative velocity error as 
ly increases. The density gradient (p(i + — p(i,j))/5 is approximately equal to the 
analytical value. No. 2 gives almost the exact density distribution. 

Of course, the achievement of second-order accuracy for Poiseuille flow does not necessarily 
mean a second-order accuracy for any flow, in next section, a 3-D duct flow will be studied. 



4.3 Stability Issue 

Another important issue is the stability related to boundary conditions. It is found that the com- 
bination of bounceback without collision on stationary walls with equilibrium distribution at flow 
boundaries gives the best behavior on stability. Once any boundary condition or flow boundary 
condition in any of the schemes in [|, §, § |(| |T| or in this paper is used, the maximum Re num- 



ber is reduced dramatically. For example, in the simulation of Poiseuille flow with bounceback 
at the wall and equilibrium scheme with velocity inlet and density outlet conditions (in the case, 
if density are prescribed at both inlet and outlet, the velocity is significantly smaller than the 
intended value for high Re flows) and with Ix = 16, ly = 8, uq = 0.1, the maximum Re is 500. 
The maximum Re reduces to 63, 56 respectively for density inlet/outlet condition No. 1, No. 2 
with bounceback boundary condition on the walls. The maximum Re further reduces to 42, 12 
respectively for density inlet/outlet condition and boundary conditions in this paper and in ||. 
It is also noted that when the parameter are close to the region of instability, the simulation may 
have unusual large errors. On this account, the half-way wall bounceback is safer. 
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5 Flow Boundary Conditions and Results for the 3-D 15-velocity 
LBGK Model 



Since 3-D model is needed in practical problems, this section will discuss the pressure or velocity 
flow boundary condition for the 3-D 15-velocity LBGK model d3ql5 and an incompressible model 
d3ql5i similar to d2q9i and present some simulation results. The model d3ql5 is based on the 
LBGK equation Eq. (||) with i = 0, 1, • • • , 14, where ej, i = 0, 1, • • • , 14 are the column vectors of 
the following matrix: 



E 



1 





-1 
1 





-1 
1 



1 
1 

-1 1 



and ej, i = 1, • • • , 6 are classified as type I, ej, i = 7, ■ ■ ■ , 14 are classified as type II. The density per 
node, p, and the macroscopic flow velocity, u = (u x ,u y ,u z ), are defined in terms of the particle 
distribution function by 



n 
i=0 



14 



P, 



5> 

i=l 



pu. 



(33) 



The equilibrium can be chosen as: 



Jo 



(eg) 



1 

-/;- -pu-u, 



(eg) 



(eg) 



-p + -pe { ■ u + -p(ei ■ u) 2 - 

11 1 , 

p + — pe-i • u + — p(ei • u) 



pu-u, i G I 



(i I 



24' 



16' 



1 

48 



pu • u. i £ II 



(34) 



The model d3ql5i is constructed from these formulas in a similar way as in d2q9i. 

The flow to be studied in the 3-D case is the square duct flow with x— direction being the flow 
direction. It is clear to give a projection of the velocities in the xz plane as shown in Fig. 1. The 
macroscopic equations of the model is the same as Eqs. (§],[|) with c 2 s = 3/8, and v = (2r — 1)5/6. 

The pressure flow boundary condition proposed in 3.1 has the following version for the model 
d3ql5: take the case of an inlet node as shown in Fig. 1, the inlet is on yz plane. After streaming, 



fi, (i = 0, 2, 3, 4, 5, 6, 8, 10, 12, 14) are known, Suppose that pi n , u y 



are specified on the 



inlet, we need to determine fi,i = 1,7,9, 11, 13 and u x from Eqs. (|33|). Similar to the derivation 
in d2q9, u x is determined by a consistency condition as: 



PinU-x — Pin [fo + h + h + h + k + 2(/ 2 + h + fw + fl2 + fu) 
The expression of x— momentum gives: 

H + h + fd + hi + /l3 = PinU x + {f 2 + fs + flO + /l2 + fu), 



(35) 



(36) 



If we use bounceback rule for the non-equilibrium part of the particle distribution fi, (i = 
1,7,9,11,13) to set 

fr = f* + i + (ft q) -ft + q h (37) 
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then Eq. (j36|) is satisfied, and all fi are defined. In order to get the correct y—, z— momenta, 
we further fix this fi (bounceback of non-equilibrium fi in the normal direction) and modify 
/7,/9,/n,/i3 as in @: 

fi ^ fi + ^eiy5y + je iz 5 z . i = 7, 9, 11,13 (38) 

This modification leaves x— momentum unchanged but adds 5 y , S z to the y—, z— momenta respec- 
tively. A suitable choice of 5 y and S z then gives the correct y— , z— momenta. Finally, we find: 



fi = fi+i + ^PinU x - ^[eiy(h - h) + e iz (f 5 - f G )], i = 7, 9, 11, 13 (39) 



There is no special treatment at the wall of the inlet/outlet if bounceback is used there. Modifi- 
cation of the flow boundary condition for d3ql5i is similar to d2q9i. 
The velocity flow boundary condition can be derived similarly. 

Simulations on 3-D square duct flow are performed on d3ql5i and on d3ql5 using the pressure 
flow boundary condition. Only the result of d3ql5i is reported, the result of d3ql5 is similar but 
the errors are greater. The analytical solution of a flow in an infinitely long rectangular duct 



-a < y < a, —b < z < b, with x being the flow direction is given by [19| 



u x (y,z)= 1 -^(-^) V (-lf-^[l - cosh(^/2a) cos(^/2a) 

Vi " ' fj/K dx i=1 g 5 ... cosh(iirb/2a) i i 3 y 1 

We use a = b = 2 in the simulations. Results of the following boundary conditions are reported: 

(1) I/O No. 1 (the I/O condition discussed above) with the half-way wall bounceback at walls; 

(2) I/O No. 2 (the I/O condition in |2|) with the half-way wall bounceback at walls; 

(3) I/O and boundary conditions in [|]. 

(4) I/O and boundary conditions in [ 10(1 - 

(5) I/O No. 2 with the original bounceback, which is the bounceback at wall nodes without 
collision (only for the middle example). 

The I/O No. 1 with solid boundary condition in [10] is also tried, the result on velocity is 
very close to that in case (4) and is not reported. Again, we fix r and the Reynolds number, 
halve 5 each time. Simulations with different Re and r < 1.3 are performed. Three examples 
are reported on Table II: (1) r = 0.6, Re = 10, (2) r = 0.8, Re = 5, (3) r = 1.1, Re = 0.2, 
and all use po = 5. The quantity Tol in the 3-D version of Eq. ( |30| ) is set as 10~ 8 . Define 
Ix, ly, Iz to represent the number of lattice steps in x— and y— , z— directions respectively, we use 
Ix = 8, 16, 32, 64, ly = Iz = 4, 8, 16, 32 respectively to do the simulation. The maximum relative 
error is defined similar to Eq. (]3l|). From Table II, we can see that the results of the half-way 
wall bounceback give an accuracy close to second-order, so do the boundary conditions in [||, jlQl . 
Besides, the errors with the half-way wall bounceback are comparable with that using the I/O 
and boundary conditions in [||, [TO]. On the other hand, the original bounceback introduces an 
error of first-order throughout the flow (LI error has a similar first-order behavior). Thus, the 
half-way wall bounceback has an essentially different behavior than the original bounceback. The 
calculated density has smaller errors than that of the velocity, and the density is less uniform 
across a section in the duct than in 2-D case. 
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It is noted that the orders of convergence in the case of 3-D duct flow with all boundary 
conditions considered are not as good as in 2-D Poiseuille flows. For example, the ratios of errors 
at the highest resolution are not very close to 4 in some cases (the order of convergence may be 
close to 2 in the cases because of an irregular large error ratio at a coarser resolution). It looks 
that in 3-D duct flow, the four edges pose additional difficulties to boundary conditions. Even 
with forcing, the density is not uniform in a cross section for the half-way wall bounceback or 
for boundary conditions in [Q, [HJ]. Nevertheless, the order of convergence for 3-D duct flow is 
still close to 2. The half-way wall bounceback has a weaker convergence when r > 1 while the 
boundary condition in [0] performs better as r > 1 but worse as r < 1. On the consideration of 
simplicity in implementation and superior stability behavior of the half-way wall bounceback, we 
think that it deserves a serious consideration in a LBGK simulation. 

We would like to point out that a second-order accuracy of the half-way wall bounceback in 
the flows considered does not imply a second-order accuracy for any flows. The statement can be 
applied to other boundary conditions as well. One is encouraged to do some tests on a simplified 
flow of the type of flows to be simulated. 

6 Discussions 

The major results in this paper are: first, flow boundary conditions can be treated in a similar 
way as wall boundary conditions, and a new way to specify flow boundary conditions based on 
bounceback of non-equilibrium part is proposed. For the test problem of Poiseuille flow with 
pressure or velocity inlet/outlet conditions, the new method recovers the analytic solution within 
machine accuracy. Second, we reconsider the the half-way wall bounceback boundary condition 
for stationary walls. It is very easy to implement. If being applied with flow boundary condition in 
this paper or in [||, it is approximately of second-order accuracy for the 2-D and 3-D channel flows 
with r being less than or close to one. The magnitude of the error is comparable with that using 
some published boundary conditions. Hence, the the half-way wall bounceback is recommended 
for stationary walls. For flow boundary conditions in simulations of small to moderate Re numbers, 
one may consider the schemes in this paper or in ||]. For the cases considered, these flow boundary 
conditions give very close results in velocity. The scheme in [Q] gives a better density distribution 
but a weaker stability behavior. One can also use the equilibrium distribution scheme if r = 1 
can be used. However, for simulations of large Re flows, one may have to use the equilibrium flow 
boundary condition with velocity inlet condition and with r close to 0.5. In that situation, the 
accuracy is only first-order. 
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Table I. maximum relative errors for the 2-D Poiseuille flow with pressure being specified at 
inlet and outlet for three cases: (1) Re=10, r = 0.6, (2) Re=10, r = 0.8, (3) Re=l, r = 1.1. The 
results of two pressure flow boundary conditions I/O Nos. 1,2 with the half-way wall bounceback 
(HWWBB) and the result of the I/O and boundary conditions in [2] are given. In each box, the 
upper figure is the error, and the lower figure is the ratio of two consecutive errors. The last 
column shows the order of convergence using the least-squares fitting. The symbol (-2) represents 
KT 2 . 







lx 


8 


16 


32 


64 


128 


order 






ly 


4 


8 


16 


32 


64 








Uq 


0.8333(-l) 


0.4167(-1) 


0.2084(-l) 


0.1042(-1) 


0.5208(-2) 








I/O No. 1 


0.6031(-1) 


0.1500(-1) 


0.3729(-2) 


0.9297(-2) 


0.2324(-3) 


2.005 


Re 


= 10 


HWWBB 


4.021 


4.023 


4.011 


4.000 






r = 


0.6 


I/O No. 2 


0.5917(-1) 


0.1479(-1) 


0.3699(-2) 


0.9265(-2) 


0.2352(-3 ) 


1.995 






HWWBB 


4.001 


3.998 


3.992 


3.939 










I/O, B.C. 






unstable 










in [2] 


















u 


0.2500 


0.1250 


0.6250(-l) 


0.3125(-1) 


0.1563(-1) 








I/O No. 1 


0.3276(-l) 


0.8319(-2) 


0.2054(-2) 


0.5111(-3) 


0.1276(-3) 


2.003 


Re 


= 10 


HWWBB 


3.938 


4.050 


4.019 


4.005 






r = 


0.8 


I/O No. 2 


0.3250(-l) 


0.8125(-1) 


0.2032(-2) 


0.5085(-3) 


0.1283(-3) 


1.997 






HWWBB 


4.000 


3.999 


3.996 


3.963 










I/O, B.C. 
in [2] 


0.1000 
4.000 


0.2500(-l) 
4.000 


0.6250(-2) 
3.999 


0.1563(-2) 
3.987 


0.3920(-3) 


1.999 






Uq 


0.5000(-l) 


0.2500(-l) 


0.1250(-1) 


0.6250(-2) 


0.3125(-2) 








I/O No. 1 


0.5550(-l) 


0.1441(-1) 


0.3617(-2) 


0.9021(-3) 


0.2249(-3) 


1.989 


Re 


= 1 


HWWBB 


4.011 


4.010 


3.984 


3.851 






r = 


1.1 


I/O No. 2 


0.5750(-l) 


0.1437(-1) 


0.3594(-2) 


0.8984(-3) 


0.2246(-3) 


2.000 






HWWBB 


4.001 


3.998 


4.000 


4.000 










I/O, B.C. 


0.5000(-l) 


0.1250(-1) 


0.3125(-2) 


0.7812(-3) 


0.1953(-3) 


2.000 






in [2] 


4.000 


4.000 


4.000 


4.000 
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Table II. maximum relative errors for the 3-D square duct flow with pressure being specified at 
inlet and outlet for three cases: (1) Re=10, r = 0.6, (2) Re=5, r = 0.8, (3) Re=0.2, r = 1.1. The 
results of two pressure flow boundary conditions I/O Nos. 1,2 with the half-way wall bounceback 
(HWWBB), the results of the I/O and boundary conditions in [2,10] are also given. The result of 
I/O in [2] with original bounceback is given for the second case. In each box, the upper figure is 
the error, and the lower figure is the ratio of two consecutive errors. The last column shows the 
order of convergence using the least-squares fitting. 





lx 

iy,iz 


8 
4 


16 

8 


32 
16 


64 
32 


order 


Re = 10 
r = 0.6 


u 


0.8333(-l) 


0.4167(-1) 


0.2083(-l) 


0.1042(-1) 




I/O No. 1 
HWWBB 


0.4028 
3.822 


0.1054 
3.844 


0.2742(-l) 
3.762 


0.7289(-2) 


1.931 


I/O No. 2 
HWWBB 


unstable 


I/O, B.C. 
in [2] 


unstable 


I/O, B.C. 
in [10] 


0.2210 
3.971 


0.5565(-l) 
4.304 


0.1293(-1) 
4.128 


0.3132(-2) 


2.053 


Re = 5 
r = 0.8 


u 


0.1250 


0.6250(-l) 


0.3125(-1) 


0.1563(-1) 




I/O No. 1 
HWWBB 


0.1382 
3.472 


0.3980(-l) 
4.059 


0.9805(-2) 
4.106 


0.2388(-2) 


1.959 


I/O No. 2 
HWWBB 


0.1371 
3.747 


0.3659(-l) 
3.959 


0.9243(-2) 
4.001 


0.2310(-2) 


1.966 


I/O, B.C. 
in [2] 


0.3397 
3.552 


0.9563(-l) 
4.517 


0.2117(-1) 
3.688 


0.5741(-2) 


1.984 


I/O, B.C. 
in [10] 


0.8567(-l) 
5.552 


0.1543(-1) 
3.427 


0.4502(-l) 
3.523 


0.1278(-2) 


1.998 


I/O No. 2 
bounceback 


0.6539 
1.955 


0.3345 
2.178 


0.1536 
1.936 


0.7935(-l) 


1.025 


Re = 0.2 
r = 1.1 


u 


0.1000(-1) 


0.5000(-2) 


0.2500(-2) 


0.1250(-2) 




I/O No. 1 
HWWBB 


0.2091 
3.199 


0.6537(-l) 
3.598 


0.1817(-1) 
3.780 


0.4807(-2) 


1.818 


I/O No. 2 
HWWBB 


0.2114 
3.279 


0.6448(-l) 
3.608 


0.1787(-1) 
3.772 


0.4737(-2) 


1.829 


I/O, B.C. 
in [2] 


0.3109 
4.017 


0.7740(-l) 
3.937 


0.1966(-1) 
4.009 


0.4904(-2) 


1.994 


I/O, B.C. 
in [10] 


0.2070 
4.162 


0.4973(-l) 
3.894 


0.1277(-1) 
3.822 


0.3341(-2) 


1.982 
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8 Figure Caption 

Fig. 1, Schematic plot of velocity directions of the 2-D (d2q9) model and projection of 3-D (d3ql5) 
model in a channel. In the 3-D model, The y— axis is pointing into the paper, so are velocity 
directions 3,7,9,12,14 (they are in parentheses if shown), while the velocity directions 4,8,10,11,13 
are pointing out. Velocity directions 3,4 have a projection at the center and are not shown in the 
figure. 
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